#include "Multigrid_Solver_1D.h"
double basic_f2rd(double * P)
{
    return  -(sin(P[0])-cos(P[0])*cos(P[0]))*exp(sin(P[0]));

};

double basic_f(double *P)
{
    return exp(sin(P[0]));
};

int main()
{
int n=128;

double _a[1]={0};
MG_1D case1;


int _Rtype;
int _Itype;
std::cout<<"Give me n and Rtype and Itype,:"<<std::endl;
std::cin>>n;
std::cin>>_Rtype;
std::cin>>_Itype;

case1.Initialize(n,basic_f2rd,basic_f);
case1.setRestircType(_Rtype);
case1.setIntpType(_Itype);
case1.setSmoothType(2);

VectorXd u_,f_,r_,r_origin,acu;
acu.setZero(n+1);
for(int i=0;i<=n;i++)
{   
    _a[0]=1.0*i/n;
    acu(i)=basic_f(_a);
};
double acu_norm=acu.lpNorm<Eigen::Infinity>();
double h;

u_=case1.get_u_origin();
f_=case1.get_f_origin();
h = case1.get_h();
r_origin=case1.ResidualSolver(u_,f_,h);
double r_origin_norm=r_origin.lpNorm<Eigen::Infinity>();
int cnt;
for (cnt=1;cnt<=100;cnt++)
{
    u_=case1.Vcycle(u_,f_,h);
    r_=case1.ResidualSolver(u_,f_,h);   
    if (r_.lpNorm<Eigen::Infinity>() <= r_origin_norm*1e-8)
    break;
};

// int cnt;
// for (cnt=1;cnt<=1;cnt++)
// {
//     u_=case1.FMGcycle(f_,h);
//     r_=case1.ResidualSolver(u_,f_,h);
//     wucha=u_-acu;
//     if (wucha.norm()/acu_norm <= 1e-3)
//     break;
// };

//  std::cout<<"The residual is\n"<<r_<<std::endl;
//  std::cout<<u_<<std::endl;

// std::cout<<"Residual=\n"<<r_<<std::endl;
std::cout<<"n= "<<n<<std::endl;
 std::cout<<"cnt= "<<cnt<<std::endl;
//  std::cout<<wucha<<std::endl;
std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;
// for(int i=1;i<2;i++){
//    case1.Vcycle(); 
// }
// solution1=case1.Vcycle();
// // solution2=case2.Jacobi();
// // acu.resize(solution1.rows());
// for(int i=1;i<n;i++)
// {
//     // acu(i-1)=exp(sin(1.0*i/n));
//      acu(i-1)=sin(M_PI*i/n);
// };
// double p=(acu).norm();
// std::cout<<solution1<<std::endl;
// std::cout<<acu<<std::endl;

// // std::cout<<"Jacobi:"<<std::endl;
// std::cout<<((solution1-acu).norm())/p<<std::endl;

};